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Abstract 

Using an operatorial formalism, we study the Kramers equation and its ap- 
plications to numerical simulations. We obtain classes of algorithms which 
may be made precise at every desired order in the time step e and with a set 
of free parameters which can be used to reduce autocorrelations. We show 
that it is possible to use a global Metropolis test to restore Detailed Balance. 
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I. INTRODUCTION 



Quantum field theories can be numerically simulated on a lattice to obtain non pertur- 
bative informations. A definite continuum limit not always exists, but for many interesting 
theories asymptotic freedom allows for the extraction of continuum physics from the lattice. 
The algorithmic problem of the Monte Carlo approach amounts to generate many statis- 
tically independent configurations of the fields. These configurations must be distributed 
according to a definite weight which depends on the specific action of the model. We have a 
great freedom in building such a procedure. A possible strategy is based on the idea of ob- 
taining the final target distribution via a diffusion process in the configuration space which 
gives asymptotically infinite configurations properly distributed. A numerical simulation 
of the above continuous diffusion process leads naturally to approximated schemes for the 
generation of the equilibrium ensemble. Indeed, solving the diffusion process is like inte- 
grating a differential equation and some kind of extrapolation in the integration step must 
be done at the end to get the exact results [|. On the other hand, exact algorithms may 
be obtained with the introduction of a Metropolis test, namely a global recall of the action 
that we are simulating. Depending on the problem under study one of the two choices may 
be preferable. In this work we show how one can generalize the simplest Langevin process 

^ In 1^ we have studied a different implementation of the diffusion process. A diffusive process 
can be realized on a discrete space-time. The correct diffusion equation in the continuum is 
obtained only at the end with a careful limit in both the spatial and temporal steps. Ref. 
showed the applicability and power of this method in the framework of the First-Exit-Mean-Time 
problems. Space becomes discrete and at every sweep the degrees of freedom are modified to 
qn+i = qn + ^a„s-a„ where {bq} identifies all the possible directions on the lattice and Aq, is the 
space step in the a-th direction. The relevant direction a is chosen with a definite space-dependent 
probability Pq, such that as t — > -|-oo with A„ ^ a definite distribution for the configurations in 
{qn} is obtained. 



2 



to more complicated and efficient schemes. We also show which global correction must be 
performed in order to restore the Detailed Balance condition. 

In Sec. |1| we recall the basic ideas of the Monte Carlo approach to numerical simula- 
tions. In Sec. |IT1| we introduce formalism and notation discussing the Langevin and the 
Horowitz algorithms. In Sec. ^ we show their natural extensions and discuss in Sec. |V] 



their autocorrelations in the continuum. In Sec. ^ we show exact versions of these schemes. 
Some comment are collected in Sec. [VII| on the relations with existing algorithms. Finally, 
Sec. [VIII] is devoted to the conclusions. In Appendix ^ we illustrate the symplectic methods 
for the integration of the equations of motion. Finally, in Appendix P and |C| some technical 
notes are collected. 



II. MONTE CARLO SIMULATIONS 

Let us consider i?" fiat space with Euclidean measure. Given the degrees of freedom {q} 
and the action S{q) (a scalar function of the q) we want to generate efficiently samples of 
statistically independent configurations distributed according to the weight 

P(g) = e-^^^'^\ (2.1) 

where j3 is the inverse temperature of the system under study. The strategy of Monte Carlo 
simulations is to start from a random initial configuration and to evolve it along a discrete 
Markov chain determined by a given transition probability 

P{q ^ q") I P{q -> q) dq" = 1 (2.2) 

which is the probability of making a transition from q to q" . We know [0] that if the above 
Markov chain is irreducible and a stationary state W{q) exists satisfying the stationarity 
condition 

I W{q)P{q ^ q) dq = Wiq"), (2.3) 
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then the process will converge to W{q). A sufficient but not necessary condition to satisfy 
Eq. ( p.3| ) is the stronger Detailed Balance condition 

W{q)P{q ^ q") = W{q")P{q" ^ q). (2.4) 

More generally, a powerful trick is to make the action depend on some auxiliary variables ^ 
which diffuse together with the q and satisfy 

S{q) S{q,0 such that J d^e^^'^'^^ ~ e^^"^). (2.5) 

Averages of functional independent of ^ are equal to those obtained with S{q). Detailed 
Balance will be encoded by a generalized Eq. ( p.3|) of the form 



Wiq,OPiq\( --q\e) = Wiq",n")Piq",Tf ^q,Ta, (2.6) 

where is the transformed of ^ under time reversal. 

Choosing the best transition probability depends on the nature of the problem under 
study and on which kind of measurements one is interested in. Moreover, for practical 
implementations, it is relevant to know the specific architecture of the machine which one 
will use. As we have explained in the introduction, a natural idea is to take advantage of 
a continuum dynamics towards equilibrium. This leads us to consider the Fokker-Planck 
diffusion equation the Langevin algorithm and its generalizations. 

III. FOKKER-PLANCK EVOLUTION 

The standard Fokker-Planck equation for the diffusion in i?" may be written 

dP 

-^-HP, (3.1) 
where the evolution kernel for the time dependent probability distribution P{q, t) is 

H = ^Il' + fflF(g) n = -z^ (3.2) 

and F{q) is the driving force which we assume derived from a potential 



F(g) = -zUSiq). (3.3) 

Eq. describes relaxation towards the static distribution 

W{q) =e-^^^''\ (3.4) 

This can be shown in the continuum evolution as follows f^. Consider the scalar product 
(we restrict ourselves to real functions) 

(/, h) = I {dq)w fiq)h{q) {dq)w = dq [Wiq)]"' , (3.5) 

then we have 

(/, Hf) = J {dq)w /(g)n jin + ^F(g)| /(g) = ~J dqW{q) {Uf{q) [Wiq)]-'}' . (3.6) 

W{q) is an eigenstate of H with zero eigenvalue. From the above equation it follows that if 
Hf = then / ~ W{q) whence W{q) is the unique stable mode. Moreover from Eq. ( |3.(j| ) we 
deduce that all the other eigenvalues of H are negative. Therefore, by expanding a generic 
solution u{q, t) of Eq. (|3TT| ) in a basis of H eigenvectors we show that u as t +oo. 

Consider now the Markov chain defined by the following real transition probability 

P{q - q) = {q"\ exp (-eH) \q') {q\K{q)\q), (3.7) 

where e is an auxiliary simulation time step. 

A necessary condition for W{q) to be a stationary distribution of this process is the 
differential equation 

K{q)W{q) = 0. (3.8) 

On the other hand, Detailed Balance may be written as the operatorial equation 

K{q)W{q) = W{q)K^{q). (3.9) 

Of course we cannot evaluate exactly the matrix elements of K. This would be equivalent to 
solving exactly the quantum mechanical problem associated to the path-integral quantum 
extension of S{q). 
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Many algorithms arise in trying to approximate the kernel K. We consider now some 
well known proposals. The basic idea is to factorize approximately the kernel K and to read 
the factorization as a succession of elementary steps in the update of the old Markov state. 



A. Langevin algorithm 

The standard Langevin algorithm splits H = Hi + H2 with 

H2 = tUF{q). 



(3.10) 



Then 



(3.11) 



The sum over ^ can be read as the composition of two successive updates realized with 
each of the two operators Hi and H2- They have trivial matrix elements in the coordinate 
representation [Q. 



27re 



n/2 



exp 



(3 II / 2 



(g"|e-^^|g')=5(g"-g(e)) 



(3.12) 
(3.13) 



where q{t) is defined by 



g(0)=g' q{t) = F{q). 



(3.14) 



The resulting explicit update is 



(3.15) 



where ^ is a Gaussian random number with zero mean and unit variance. 

This procedure introduces a systematic error in the construction of the statistical sample 
since the static distribution is now the zero mode of the operator H defined by 
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exp {—eHi) exp {—eH2) = exp (^—eH^ 



H = H + e5iH + 



(3.16) 



Solving perturbatively in e the equation 



(H + eSiH + ...){W + eSiW + ...) = 



we find a non trivial 6iW which affects to order 0(e) all the averages on the sample 



(3.17) 



W + e5iW = exp -/3 



s + -(v^s-(3{vsy 



(3.18) 



Therefore, Detailed Balance holds exactly only with respect to H. Of course, we are not 
able to build explicitly the function q{t). At this order we can use 



q =q + eF{q) + Jj ^ 



(3.19) 



and the only consequence is a different explicit form of the correction term which becomes 
now 



W + eSiW = exp <^ -p 



s+'-h's-^ivsr 



(3.20) 



1. Smart Monte Carlo 



In this section we show how the structure of the Detailed Balance violations in the 
Langevin algorithm (which the Smart Monte Carlo method tries to correct) suggests 
the extension of the dynamics in an enlarged {q,p) space. In this wider space the diffusive 
process modifies some approximation of a canonical Hamiltonian evolution. The Langevin 
update is (e — > e^/2) 

^2 



q =q+-F{q) + ei 



(3.21) 



corresponding to the transition probability in Eq. (|3.12| ). The interesting point is that the 
parameter controlling Detailed balance violations 

W{q){q\K{q)\q) 



X 



W{q){q\K{q)\q) 
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(3.22) 



may be written 

X = exp {H{q",p") - H{q , (3.23) 
where H{q,p) = S{q) + and 

P = ^, 

q" =q' + jFiq) + ep, (3.24) 

p =p +-[F{q)+F{q )) , 

and this set of equations can be recast in the form of a standard leap-frog symplectic inte- 
grator (see Appendix 0) 

P = -C, 

P =P + |^(?')' 

q" = q + ep, (3.25) 
P =P+^F{q ). 

The conclusion is that the equality x = 1 can be enforced improving the approximation of 
the evolution determined by H. 



B. Hyperbolic algorithm 

This is a clever proposal which improves the naive Langevin algorithm using the 
dynamics of the Kramers equation We double the degrees of freedom q — > {q,p) and 

introduce the operator 

H = ^Ul + ^^p(-7p) + tIl,F{q) + tU,p, (3.26) 

where 7 is an arbitrary parameter. 

The proof of convergence to equilibrium given for the Fokker-Planck evolution cannot 
be applied. However, following P|, H determines a Hamiltonian flow in {q,p) space with a 
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diffusion term in the p variables. Hence if a solution of m = Hu is positive at t = 0, it remains 
positive at any time later as can also be seen from the representation of the time evolution 
of H in terms of a path-integral with positive measure. On the other hand, integrating by 
parts any eigenvector u„ with non zero eigenvalue we get 

J dq dpun = (3.27) 

and assuming some eigenvalue with positive real part we contradict positivity of the evolution 
determined by H. 

We introduce an 0(e) error by splitting H = Hi + H2 with 

Hi = ^Ul + tUp{-fp) + iUpF{q), (3.28) 

H2 = iUqP 



we can write down an approximated update for H obtained from Hi and H 
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p" = {l-e^)p +eF{q) + ^'^e, (3.29) 
q = q + ep . (3.30) 

The equilibrium distribution is 

W{q,p) = exp + S{q) + eSi{q,p,^) + ■■■)}• (3-31) 

The parameter 7 can be tuned to minimize the autocorrelation time of the observable which 
we are going to measure. Unlike the Langevin case, it is impossible to give an expression of 
Si{q,P,l) in the form of a polynomial in S{q) and its derivatives. 

On the basis of the above general arguments we should have introduced a systematic 
error of order 0(e) due to S'i(g,p, 7). This holds true for averages of generic functions 
Q{q,p) which depend on the auxiliary variables. However, if we are interested in functions 
f2(g) of q only, it is shown in that the particular form of >S'i(g,p, 7) reduces the error to 
O(e^). 

Assuming the q even under time reversal, the p variables must be are odd under T. This 
implies that the Detailed Balance condition must be written 



K{q,p)Wiq,p) = W{q,p)K^iq, -p) (3.32) 

and in our case it is violated by 0(e) terms because of Si{q,p,'-y). 

An interesting variant of this procedure is called Guided Hybrid Monte Carlo 0, this is 
an algorithm which admits an exact extension and we shall discuss it later as a particular 
case of our methods. 



IV. IMPROVED ALGORITHM 

In this section we shall describe our approach which allows us for an improvement of 
the precision up to any desired order and which is a natural generalization of the previous 
schemes. 

The key point of our method is a different operator splitting of H. Indeed, we can 
separate H = Hi„. + Hrev where 

i/,,, = ^n; + mp(-7P), (4.1) 

Hrev = ^npF(g) + tU,p. (4.2) 
The important point is that both Hirr and Hrev annihilate the static distribution 

Wiq,p) = exp + S{q)^y (4.3) 

After all, this is why the Kramers equation works in the continuum limit. Writing 

exp {-eHirr) exp {-eHrev) = exp (-e {Hirr + Hrev + eX)) , (4.4) 

the operator X is non trivial but has the useful property XW{q,p) = 0. This observation 
has the important consequence that we can consistently set to zero all the corrections in 
powers of e to W{q,p) and treat independently the irreversible and the reversible evolutions. 

We expect Detailed Balance to hold exactly: indeed, as we have remarked above, the 
two kernels 

exp {-eHirr) exp {-eHrev) (4.5) 
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leave invariant W which is therefore a steady state for the evolution determined by any 
string K of the form 

= n {-O'i^Hirr) exp {-hieHrev) , (4.6) 

i 

where and bi are constants. 

Indeed, we can check that separately 

Hirr{q,p)W{q,p) = W{q, -p)HUq, -p), (4.7) 
Hrev{q,p)W{q,p) = W{q, -p)HlM. -p), (4-8) 

(4.9) 

we conclude that the Detailed Balance condition is exactly satisfied by a symmetric splitting 
like for instance 

K = exp (^—eHirr^ exp {-eHrev) exp (^—eHirr^ . (4.10) 

Of course, we must specify how to evaluate the matrix elements of the two kernels, but this 
is not a problem. The irreversible part can be solved exactly. Writing 

= ^ K (4-11) 

02 = ^^J,(-7p), (4.12) 

then we obtain 

[ni,fl2] = -2jfli (4.13) 
and by the Campbell-Baker-Hausdoff formula 

exp {-6(^1 + VI2)) = exp {-h{e)Vli) exp (-6^2) ■ (4.14) 
Differentiating with respect to e with the condition h{0) = we get 

Me) = ^ (1 - e-^^^) . (4.15) 
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The matrix element (p", q)irr corresponds exactly to the discrete update 



p — e ^ p + 




(4.16) 



The reversible matrix element (p", q'\p , q)rev must be some approximate integration of the 



deterministic Cauchy problem 



(4.17) 



corresponding to the update 



q" ^Q{q',p',e) p" ^ P(q' ,p' ,e), 



(4.18) 



where the maps Q and P can be, for instance, an higher order symplectic integrator. We 
stress that at this level of the discussion we do not need using a canonical integrator con- 
serving the Poincare invariants, in particular the phase-space measure. The importance 
of symplectic integrators will be clear in the discussion on the exact extensions of these 
algorithms. 

We remark that the systematic error we are introducing is just the error made in approx- 
imating this reversible step. Its improvement is a well known problem which can be solved 
up to an arbitrary precision with increasing computational cost. 

As in the hyperbolic algorithm, the free parameter 7 can be tuned minimizing autocor- 
relations times, whle the integration step e must be kept small enough to control systematic 
errors. 



Following the above ideas, we can introduce schemes in which other free parameters are 
present. They can be tuned to obtain even smaller autocorrelations times. The simplest 
example is 



A. Extended schemes 



roti 



(4.19) 
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where 

H,„^lul + iUy{-jy), (4.20) 

Hrev = illpFiq) + iUgp, (4.21) 

Hrot^ia{yTip-pTiy)- (4.22) 

The operator Hrot is just a rotation in the {y,p) plane and its exact matrix elements 

{p"y"\exp {-eHrot) \py') (4.23) 

corresponds therefore to the update 

p" — p cos(Q;e) + y' sin(Q;e), (4.24) 
y = y cos(Q;e) — p sin(ae). 

A symmetric kernel like 

exp {-eHirr) exp {-eHrot) exp {-eHrev) exp {-eHrot) exp {-eHirr) (4-25) 

satisfies exactly the Detailed Balance condition with respect to the equilibrium distribution 

W{q,p, y) = e-'^(^(p'+3/')+5(9)). (4.26) 

The parameters 7 and a must be tuned, but wc arc unable to provide a general rule for 
their choice. However we shall show that the additional y variables improve the optimal 
autocorrelation of the algorithm. 

An interesting feature of this scheme is that if we let y reach equilibrium, inserting one 
exp {—eHirr) factor near exp {—eHrev) and setting 7 00, then we obtain 

y' = ^«. (4.27) 

p = pcos(ae) H — ^^sin(Q;e) 
VP 

which is just the mixing which we find in the stochastic step of the Guided Hybrid algorithm. 

It is clear that, at least formally, we reduce to this update starting from our basic {q,p) 
scheme and relating 7 to o; writing 
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exp(— 7e) = cos(Q;e) < cte < — (4.28) 

(the case in which the cosine is negative is handled easily by a slight generalization introduced 
in the following section). 

A very natural extension is obtained building in the space {q,p,yi, • ■ ■ jUn) a chain of 
rotation operators acting in the space of the auxiliary variables 

Hrev = i^pF{q) + iUgp, (4.29) 
Hrot,l = ^Ol{ylIlp-pUyJ, (4.30) 

(4.31) 

Hrot,N = iOn {yN^yN--L " VN-J^yN) > (4-32) 

Hirr = ^ + i^vA-lVN), (4.33) 
where the angles 9i are free. 



B. Some remarks 

1. Arbitrary rotations 

The reason why it is so natural to introduce the rotation operators H^^t is the following: 
if p be a A'"-component vector, then the discrete update 

p" = xCp + Vl^i C e 0(N, F) (4.34) 

generates a normally distributed succession {p}: Pip) — exp 

The connected component of this update may be written in terms of the exponential of 
e times a differential operator by introducing all the generators of so(N, IR) 

Hrot = i ©ii Pi rip. Qij = -Qji. (4.35) 

The extended scheme of the previous section is just a particular choice of C realized as a 
chain of elementary rotations. 
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In the steps where the matrix element of the transition kernel are computed exactly, 
the parameter e does not control an approximation any more. Indeed, our convergence 
analysis deals with the limit e ^ irrespectively of 76, Qije. On a formal level, one can take 
parameters 7 or Qij not 0(1) with respect to e. From a practical point of view this means 
that we can choose C independent of e; intuitively when the dynamics is switched off one is 
left with a free rotation in the p space0. 

The vector p may be the full set of generalized momenta in a single site. However, it can 
also mix momenta from different sites. This introduces long distance reshufflings and allows 
for pair-exchanges which are represented by orthogonal matrices. 

2. Random parameters 



In the Hybrid Monte Carlo [11,0], better autocorrelations may be obtained by random- 



izing the trajectory length. Here, a probability distribution for parameters like 7 may be 



^ An interesting example of this mechanism if the following. We could try looking for an extension 
of the basic update 

/ = ap + /?C (4.36) 

to the case a S C and (3 to be determined. However, by complexifying p and representing them 
as 2-vectors 

// 

p = u + iv ^ = rj + ip i (4.37) 

we have 

p' = \a\C'p + I3i C = ^ (Reo - icj2lma) G 0(2, F) (72 is a Pauli matrix (4.38) 

\a\ 

and choosing (3 = yjl — \d^\ we fall in the above case, since C may have a limit different from 1 as 
lal — > 1. 
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introduced. However, from preliminar numerical investigations we have indications that this 
device should not be a major improvement. 



3. Dynamical Fermions 



Consider a model with djTiamical fermions. As is well known, in such a situation, the 
most expensive part of a lattice simulation is the inversion of the fermionic propagator. The 
time spent during this step becomes the natural unit of time. A fundamental parameter is 
then the number of inversions per sweep which we shall call TZ. A common implementation 
of dynamical fermions [|ll|] replaces the fermionic degrees of freedom with bosonic fields x 
which are updated by a heat-bath step. The refreshment frequency is relevant here. As 
an example, let us integrate the equations of motion by the simplest leap-frog scheme of 
Eq. ( [A18| ). If we refresh the x fields every k sweeps and if N^d is the number of molecular 
dynamics steps in each sweep, then 



In the Hybrid Monte Carlo one usually takes N^d ^ 1 and k = 1. In the schemes we have 
discussed one has always N^d = 1 so that taking A; = 1 is twice slower than a large k. 

V. DISCUSSION OF CORRELATIONS 

Consider a function Q{q) of the dynamical variables. The evolution in time of the sample 
average of Q is given by 



Consider the continuum limit in which e ^ and let us compare the evolution of (g) in the 
free theory where 



n = N^d + T- 



(4.39) 




(5.1) 




(5.2) 
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for the three cases discussed above. We restrict the discussion to the case in which all the 
tunable parameters are independent of e (see the discussion in paragraph ( |1V B We have 
explicitly 

-Hlyp = ^dl- -fpdp - u^qdp + pdq, (5.3) 

1 ^ 
~^Lp,y) = ^ - IVNdy^ + J2 {Vidy.-i - Vi-idy,) - uj\dp + pdq, 

iVo =P)- 

The eigenvalues of the time evolution of (g) give directly the intrinsic autocorrelation of the 
Markov chain in the e — >^ limit. 

For the Langevin algorithm, we have a decay with eigenvalue A = —u^ and no possibility 
of tuning any parameter. 

In the hyperbolic case, we have two eigenvalues X± and 

minmax(A+, A_) = —u (5.4) 

7 

obtained at 7 = 2lj. Dimensions are different with respect to Langevin, because the Kramers 
equation is second order; afterall, this is why the hyperbolic algorithm is better than the 
Langevin one. 

In the (g, p, y) extended scheme the roots of the secular equation Aj are tedious functions 
of 6i and 7. We remark that a possible choice is to choose 7 and 6i in such a way to have 
the maximum possible degeneracy in A. Some explicit results for low values of are shown 
in Tab. | where all the variables has been made adimensional dividing out by u. One can 
improve the autocorrelation of (q) of at least a factor ^/TT in the = 5 case. It seems 
reasonable that this mechanism can be improved with no limit introducing more and more 
auxiliary variables (up to < 6 Tab. H satisfies Aat = —\^2N + 1) and we must balance the 
increasing computational work needed to tune all the associated new parameters and the 
gain in autocorrelation of the sample. 
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If we take into account the details of the discrete dynamics at e > we obtain hnear 
recursive equations giving the evolution of the momenta instead of differential equations. 
Their characteristic polynomials determine the e dependent relaxation eigenvalues. For e 
small enough the qualitative picture which one finds in the continuum does not change. 

Of course, we cannot be too optimistic. In a realistic model there will be many relevant 
frequencies. An optimal tuning for all of them will average on the their separate behaviours. 
Since we cannot give analytical estimates for autocorrelations in interacting models, explicit 
numerical tests are needed. 



A. Curved Space 

We are mainly interested in the generalization from flat space to the case of the curved 
manifold of a Lie group G. We must endow the configuration space with a definite metric 
and write down the corresponding covariant Fokker-Planck or Kramers equation. The most 
convenient and natural choice is the left and right invariant metric ds^ = Tr((if/ dW). 
Its Laplace-Beltrami inavriant operator is built by substituting derivatives with left- 
derivatives on the group defined by 

/(exp(ze"T")^7) = f{g) + e^VJ{g) + ©(e^). (5.5) 

This choice is convenient because da was the generator of translations in flat space whereas 
on the group is simply the generator of translations by left multiplication 0. In the 
Langevin case, the diffusive term is the heat kernel 

(^?1e-^^^|/) = (l|e-^^^|(^')"'/), (5.6) 
namely the matrix elements of the Laplace-Beltrami operator on G. For SU{N) a compact 



expression can be given flQ] which however is gaussian in the Lie parameters only if g' is 



^We recall that dealing with we can integrate by parts with respect to the Haar measure 
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near to g (the weak coupling limit of QCD). The advantage of the Kramers equation is 
that it is associated to a diffusion on i?" x G where n = dimC, the p variables being fiat. 
The molecular dynamics on the group is more complicated than in the flat case, but no 
curved heat kernel is needed here because the Hamiltonian evolution on G is modified by a 
diffusion in its flat Lie algebra. 



Until now we have considered methods which are not exact and need an extrapolation to 
obtain the exact e ^ hmit. Moreover, finite e effects can be studied by the analysis of the 
effects of the correction terms in the equilibrium action and it is not always clear if critical 
behaviour is unchanged if this terms are present. Nevertheless, they have some advantages 
over exact methods; on a formal level, they can deal with complex actions and they do not 
suffer the decrease of efficiency with increasing volume of the Hybrid Monte Carlo algorithm 
due to the acceptance rate. 

Anyhow, adding a global accept/reject test, we can write exact extensions which satisfy 
the Detailed Balance condition. Let us consider the minimal scheme of the improved 
Hyperbolic algorithm; let us introduce an irreversible transition probability H, a reversible 
one T and an accept/reject test associated with the acceptance probability A. We want to 
prove Detailed Balance in full generality so we point out which properties of these compo- 
nents we need. 

The irreversible transition probability must satisfy 



VI. EXACT EXTENSIONS 



H(p ) 



H(-p' 




(6.1) 




(6.2) 



VF(g,p )n(p — >• p ) 



W{q,p)Ii{p' ^p), 



(6.3) 



where we have put /3 — 1 and W{q,p) = exp{—S{q,p)) is the equilibrium state. 



The properties of the reversible transition probability are 
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J T{q\p' -> q",p") dq" dp" = 1, 

m/ ' ' II II . , II II I 1^ 

T{q ,p ^ q ,p ) =T{q ,-p q , -p 



The requirement on A is the following: define 



K{q',p' q",p") = A{q,p q",p")T{q,p' q" ,p 



I I II II . „ , I I II II 



then we want 



/ I ^ ^ ^ y f I II II 



II II ^ ^ ^ , II II I / N 



W{q ,p)K{q ,p -> q ,p ) = W{q ,p )K{q , -p q , -p 

An explicit choice of 11, T and A which satisfies all the above conditions is 

1 



Il{p p ) = c exp 



-{p —xpY 



2(1 

T{q',p' q",p") = 6{q" - Q{q\p'))6{p' - P{q,p)), 



A{q',p' q",p") = min{l, W{q' ,p')/W{q ,p)}, 



(6.4) 
(6.5) 



(6.6) 



(6.7) 



(6.8) 



where x is any real number and Q, P are the maps associated to a symplectic integrator. 
The integrator must be canonical to avoid any Jacobian in verifying Eq. (|6.5|). 



A. {q,p) scheme 

Our extension of the (g, p) scheme is the following 

{q\p') ^ {q\pi), 
{(l\pi) ^ (?1,P2), 

accept : (g",p3) = (^1,^2) with P = A{q',pi ^1,^2), 
reject : (^",^3) = (-pi, g') with P = 1 - A{q' ,pi ^ qi,p2), 

(g",P3) ^ {q",p")- 



(6.9) 



The full transition probability of the above algorithm is 
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P{q ,p q ,p ) = 

J n(p' pi)K{q',pi q'\p2)Il{p2 p") dpi dp2 + 

+ J n(p' pi){l - A{q',pi qi,P2)}T{q' ,pi -> ^1,^2) x (6.10) 

xn(-pi p")6{q' - q") dpi dp2 dqi. 

Using the properties of 11, T and A we verify the normahzation condition 

P{q ,p ^ q ,p ) dq dp =1. (6-11) 

In Appendix |B] we show that 

W{q\p')P{q' ,p' q'\p") = W{q" ,p")P{q'\ -p" g', -p), (6.12) 
hence Detailed Balance is demonstrated. 

B. {q,p,y) scheme 

We restrict ourselves for simplicity to the = 1 case. If we assume q even under time 
reversal, it follows that also y is even whilst p is odd. We expect to find a Detailed Balance 
equation of the form 

. I t I . I I t II ir II . , II II II . , II II II I II . . 

W{q ,p ,y )P{q ,p ,y ^ q ,p ,y ) =W{q ,p ,y )P{q ,-p ,y ^q,-p,y). (6.13) 

We must remember that the matrix element of Hrot in the (p, y) plane is exact and leaves 
invariant the sum of squares in the action. Hence we have 

W{q ,p' ,y')Trot{q\p ,y' ^q",p'\y") = W{q',p",y")Trotiq\p\y' ^q",p",y"). (6.14) 
The scheme we propose is 

{(i\p\yi) ^* iQ\puy2), 
i(i',pi,y2) {quP2,y2), 
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accept : {q\pz) = (^1,^2) with P = A{q',pi ^1,^2), (6.15) 
reject : (g",P3) = {-Pi,q') with P = 1 - A{q',pi 51,^2), 

/ " \ Trot f " " N 

{q ,P3,y2) {q ,p ,^3), 

II II X n / " " "n 

[Q ,p ,1/3) (g ,p ,1/ )• 

In Appendix y we show that even for this scheme, the Detailed Balance condition in 
Eq. ( imi) holds. 

VII. COMPARISON WITH HYBRID MONTE CARLO 

Let us first consider the (g, p) scheme. This algorithm has been tested on compact QED 
in where its performance has been compared with the global Hybrid Monte Carlo. 

In the 7 ^ 00 limit, we recover exactly the Hybrid Monte Carlo algorithm with a single 
molecular dynamics step. The quantity (7e)~^ plays qualitatively the role of Nmd and a 
finite 7 tries to reproduce a N^d > 1 dynamics by taking into account the old momenta of 
the microcanonical evolution. 

We remark that, after integrating out the fermions, the non local nature of the action 
prevents one from taking advantages from local Hybrid Monte Carlo algorithms |T^. The 



expected scaling behaviour is discussed as follows. Assume our lattice model to satisfy L 3> 
^, where L is the lattice size and ^ the correlation length, then the acceptance probability of 
the HMC behaves as Pace ~ erfc(cA^mde^V^)- Optimal tuning is expected to give N^d ~ 1/e 
with a sweep-sweep correlation r ~ 1/e. To keep the acceptance rate constant while varying 
the volume needs the scaling e ~ V~^^'^. The Kramers algorithm has N^d = 1 independently 
on the tuning parameters. We assume that the scaling relation is modified to e ~ V~^^^; 
this assumption amounts to neglect the effects of momenta negation on the acceptance 
probability. If the optimal autocorrelations stay ~ 1/e then under our assumption, the 
HMC is expected to behave badly on larger volumes. Moreover, from the point of view of 
numerical precision, the opportunity of having N^d = 1 is welcome; indeed, in QCD, as 
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the volume (L^) grows or the quark mass (m) is decreased we must have N^d ~ L/rrv'/^ 
and some protection seems necessary to protect irreversibihty against the accumulation of 
numerical errors. 

The same qualitative discussion applies to the extended schemes. In a realistic model 
many frequencies are present and the A values of Tab. | are too optimistic. However, we 
expect an improvement of efficiency with an appropriate tuning of 6i and the same qualitative 
scaling with volume. 

VIII. CONCLUSIONS 

In summary, we have analyzed the Kramers equation approach to lattice simulations. Our 
operatorial approach is an unified framework in which the convergence analysis is straightfor- 
ward. Generalized algorithms come out easily with an increasing number of free parameters. 
They eventually can be tuned in order to minimize autocorrelations. Explicit numerical tests 
are needed. All these algorithms are particularly safe from the point of view of numerical 
precision and they are indicated when a great deal of computational effort is needed. Sub- 
sequent papers will be devoted to a numerical confirmation of the convergence and scaling 
properties of these proposals. In particular, we are interested in models with dynamical 
fermions where the numerical precision problem may be important. 
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APPENDIX A: SYMPLECTIC METHODS 

In this appendix we address the problem of the numerical canonical integration of the 
Hamilton equations. Many studies of this topic can be found in literature. Because of 
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its importance these works appear in many different research fields. We quote [p^-[T6 



as examples of papers belonging to the lattice field theory literature. However, similar 



investigations could be find for instance in [1171 ^ different context. An useful review for 



the reader may be |T^ where specific features of the problem are discussed including (i) non 
separable Hamiltonians (ii) situations in which the computation of the derivatives of S{q) 
is not expensive. 

We want to review shortly the techniques needed to build efficiently higher order ap- 
proximation schemes. Let us consider a separable Hamiltonian H = T{p) + V{q) with 
degrees of freedom. The symplectic 2-form uj"^ = dp A dq and its exterior powers {uj'^Y with 
k = 1 ■ ■ ■ N define the Poincare invariant integrals which give invariant quantities when in- 
tegrated; Liouville theorem is obtained with k = N. We want to define a numerical map of 
initial conditions {qo,Po) such that u"^ is exactly preserved and the flow of H is approximated 
in some senseQ. 

Let (g(e),p(e)) the exact evolution of {qo,Po) under H. A symplectic integrator of order 
n is defined to be a canonical transformation $ 

(^o,Po) ^ (g,P,e) (Al) 



such that 



{q{e),p{e)) (go,Po) = (go,Po) + 0(e"). (A2) 



^ We remark that canonical integrators are intrinsically more stable than the non canonical ones, 
but they do not solve completely the problem of numerical stability. One would be tempted to 
apply KAM theory to conclude that a conserved new hamiltonian exists in the numerical evolution 
of the integrator. However, KAM theory does not apply here; as we shall see we are dealing with 
a Hamiltonian with discontinuous changes in time. Indeed, in iV > 2 degrees of freedom a large 
numerical diffusion may be observed (Arnold's diffusion) even with canonical integrators. One 
usually expects this effect to be small, but in general there is no hope to catch the long time 
dynamics of H. 
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Finding symplectic integrators is easy because of the following observation. Consider a 
canonical transformation — > {qo,Po) such that {qo,Po) — > {qo,Po) as e — and that the 
new Hamiltonian is 

oo 

H{%,Pq) = hk{qo,Poy (A3) 

k=n 

then the Hamiltonian evolution of {qo,po) — > {qo,Po) is zero up to order ©(e""*"^) and Eq. ( |A2[ ) 
holds. Consider now the chain of canonical transformations 

{qi,pi) ^ {qi-i,Pi-i) ■ ■ ■ ^ (go,po), (A4) 

where the generating functions are 

Ki{qi_i,pi,e) = -qi-iPi - [aiT{pi) + biV{qi_i)] . (A5) 

The explicit transformations are 

qi = -dp^Ki = qi_i + eaidpT{pi), 
Pi-i = -dq^_^Ki =pi + ebidp^_^V{qi-i), (A6) 
Hi_i{qi_i,pi_i) = Hi + d,Ki. 



Imposing the validity of Eq. (|A3|) we obtain the desired symplectic integrators. We remark 
that the number of steps / grows faster than the order n like in standard Runge-Kutta 
schemes and finding higher order schemes by brute force is not possible. For the actual 
construction of symplectic integrators a different approach, namely the operatorial one, is 
more practical. Introducing the notation 



DgF = {F,G} = J2(^ 



dF^dG_d£dG\ 
dqi dpi dpi dqi ) 



(AT) 



the Hamilton equations are 

z = Dhz z(e) =$exact(e) ^2(0), (A8) 

where 
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$exact(e) =exp[e(A + B)] A = Dt(p) B = Dy^q) (A9) 

and we have expressed the Hamiltonian evolution in terms of the exponential map of the 
sum of two differential (non commutative) operators. The explicit expression of $exact is 
not known in general. Again, an integrator of order n is a canonical approximation $n(e) of 
exact exact to order 0{e^). Recalling Eq. ([A5|) we choose $ in the factorized form 



k 

^{e) = f{e"'^^ e"^^^. (AlO) 

i=l 

The associated map z(0) — > z{e) corresponds as before to the chain of infinitesimal canonical 
transformations 

Qi = qi-i + eCidp^^^T{pi_i), (All) 
Pi = pi-i - edidq^_^V{qi^i). (A12) 

The unknown constants {q} and {di\ may of course be determined by brute force. However, 
a key observation which simplifies computation is that for any two operators A and B we 
have 

expeA expe5 exp = exp |e(2A + 5) + [S, A]] - [A, [A, S]]) + 0(e^)| . (A13) 

If $ satisfies the time reversibility equation 

$(-6) =$-1(6), (A14) 

then $ does not have even terms 

<l>(e) = exp I e^^fcj with X2 = X4 = • • • = 0. (A15) 

Using this property we can proof that given a symplectic integrator of order 2n satisfying 
Eq. ( |A14| ), the integrator 



^2n+2{(^) = <^2n{Zie)<^2n{Zoe)<^2n{Zie) (Al6) 

is symplectic, satisfies Eq. ( |A14|) and is of order 2{n + 1) if and only if 
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2l/(2n+l) ]_ 

Since $2 exists and is obtained using 

k = 2 ci = C2 = ^ di = l (A18) 

in Eq. ( [A10| ), we can show inductively that a reversible symplectic integrator of arbitrarily 
high order 2n exists which involves 3""^ $2 factors, namely 3"~^ + 1 elementary steps. 

We remark that when high order integrators are needed, more direct approaches can save 
time. As an example, in [|l7l we find the proposal 



$M = $2(6^™) ■ ■ -$2(6^0) ■ --M^^m). (A19) 

By a straightforward computation one finds that m = 3 and m = 7 are enough for the 6th 
and 8th order integrator. The constants {wi} are given as numerical solutions of a set of 
algebraic equations. The reduced number of steps required is 8 and 16 in this case, whereas 
the scheme of Eq. ( |A16| ) needed 10 and 28 steps respectively. 



APPENDIX B: Proof of Eq. (|6J^ 

We write 



Wiq,p')P{q',p ^ q",p") = X, + X2 + (Bl) 



with 



Xi = W{q',p') J U{p' pi)K{q',pi q",p2)Il{p2 p") dpi dp2, (B2) 



X2 = W{q',p')Siq' - q")U^'\p' ^ -/), (B3) 



Xs = -Wiq',p') jYi{p -^p^)K{q ,Pi (I11P2) X (B4) 



xn(— pi — p )6{q — q ) dpi dp2 dq 



1) 



27 



where 

n(2)(p, ^ P2) = / n(pi ^ p)n(p ^ ^2) ^^p- (B5) 

Rewriting Xi as 

Xi = W{q\p') j n(/ -> p2)K{q", -p2 q, -pi)Il{pi p) dpi dp2 (B6) 
it is clear that 

,p g ,p ) = , -p ^q,-p). (B7) 

The same propery holds for X2 since 11^^^ satisfies Detailed Balance like 11 does. Finally 



,11 II lU 

- J dqi dpi dq2 dp2 6{q2 - q')6{q2 - g")n(pi p')n(pi -> -/) X (B8) 
xK{qi, -p2 q2, -pi)W{qi,p2) 

from which it is evident that 

^3(9', q",p") = Xsiq", -p" q, -p). (B9) 



APPENDIX C: Proof of Eq. ( IsTTsD 

We now proof Eq. ( |6.13| ) . We write as before 



W{q\p\y')P{q\p\y' -^q\p\y")=Y, + Y2 + Y3 (CI) 
then 

Yi = W{q',p\y') J n{y' -> yi)Trotiyi,p' y2,Pi)K{q,pi q",P2) x 

xTrotiy2,P2 ?/3,/)n(i/3 ^ v" ) dpi dp2 dyi dy2 dy^ (C2) 

and we see that 
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, / / / It ff ft. -xr / " It tl t II. / /-ie-t\ 

Yi{Q:P:y ,p ,y ) ,y ^q^-p^v)- (C3) 

Regarding Y2 

Y2 = W{q\p\y')6{q' - q") j Il{y yi)Trot{yi,P ^ y2:Pi)T{q\pi ^1,^2) x 

xTrotiy2, -Pi 2/3,/)n(|/3 y") dqi dpi dp2 dyi dy2 dy^, (C4) 



but 



j Trotiyup y2,pi)Trotiy2, -pi y3,p") dy2 dpi = 

= / {ys, -p"\e'"^°*\y2,pi){y2,Pi\e-'"^''*\yi,p) dy2 dpi = (C5) 



(1/3, -P \yi:P) = S{yi-y3)6{p +p 



whence 



Y2 = W{q',p',y')5{q' - q")5{p' +p")n^'\y' - y") (C6) 



and again 



, / / / II II II, .11 II II I I /, I r-M-7\ 

Y2\q,p,y :P :y ) = >2(? , , ^ ^q,-p,y)- (C7) 



Finally, 



>3 = W{q ,p ,y)5{q - q") j U{y' yi)Trot{yi,p' ^ y2,Pi)K{q' ,pi ^1,^2) X 

xTrotiy2, -Pi ^ 2/3,/)n(2/3 y") dqi dpi dp2 dyi dy2 dy^ = (C8) 
= 5{q' - q") j n(yi y')U{ys y")Trot{yi,P ^ y2,Pi)Trot{y2, -Pi ^ y3,p") X 
K{qi, -P2 q, -Pi)W{q' ,pi,y2) dqi dpi dp2 dyi dy2 dy^. 

Since in the motion determined by Hrot we can view y as a. coordinate it follows that 

Trot{y\p' y",p") = Trotiy", -p" ^ y\ -p) (C9) 



We can verify that changing p and p' into — p", —p and exchanging the names of yi and 
7/3 the expression for F3 does not change. 
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TABLES 







TABLE I. 


Some admissible 7, 


0j for the {q,p,y) scheme 






N 


A 


7 


^1 




03 






1 
2 
3 

4 

5 


-V3 
-V5 
-V7 
-3 

-\/Ti 


3^3 
18 

7\/TT 


2V2 
2 

4/V5 


5 

/243 
V 35 

11 

\21 


V56 

/96 
V 5 


v/33 


iVTi 
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